library(reshape2)
library(ggplot2)
library(grid)
library(gridExtra)
library(mediation)
library(stargazer)
library(psy)
library(survey)
library(lmtest)
library(Hmisc)
library(Amelia)
library(AER)
library(DirectEffects)
library(TOSTER)
options(mc.cores = 4)

# Please change this working directory to the location of the replication data
setwd("~/Replication_Data/")

# ggplot theme
source("ggplot_theme.R")

D <- read.csv("ParisData.csv", stringsAsFactors = FALSE)



##########
########## CLEAN & REORDER DATA, AND CREATE INDEXES
##########

# Respondents who have complete socio-demographics
D$allSocio <- as.numeric(!is.na(D$gender) &
                         !is.na(D$ageGroup) &
                         !is.na(D$education) &
                         !is.na(D$motherTongue) &
                         !is.na(D$region) &
                         !is.na(D$selfPlacement))

# Reorder factor levels for education
D$education <- factor(D$education, levels = c("High school or below", "College", "University degree"))

# Creating the indexes

# Create sympathy index / standardize
D$indifference <- abs(D$indifference - 11) # Reverse direction of indifference scale
D$sympathyIndex <- ((D$sympathy + D$indifference + D$compassion + D$sadness) - 4) * (10/36) # Re-scale to 0-10
D$sympathyIndex <- (D$sympathyIndex - weighted.mean(D$sympathyIndex[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$sympathyIndex[D$period == 0], D$weight.lang[D$period == 0]))

# Create anxiety index / standardize
D$anxietyIndex <- ((D$anxiety + D$upset + D$worry + D$fear) - 4) * (10/36) # Re-scale to 0-10
D$anxietyIndex <- (D$anxietyIndex - weighted.mean(D$anxietyIndex[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$anxietyIndex[D$period == 0], D$weight.lang[D$period == 0]))

D$anxietyWithAnger <- ((D$anxiety + D$upset + D$worry + D$fear + D$anger) - 5) * (10/45) # Re-scale to 0-10
D$anxietyWithAnger <- (D$anxietyWithAnger - weighted.mean(D$anxietyWithAnger[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$anxietyWithAnger[D$period == 0], D$weight.lang[D$period == 0]))

# Standardized anger indicator
D$angerStd <- (D$anger - weighted.mean(D$anger[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$anger[D$period == 0], D$weight.lang[D$period == 0]))

# Create economic threat index / standardize
D$economy3 <- abs(D$economy3 - 7) # Reverse scale direction
D$economyIndex <- (apply(D[, c("economy1", "economy2", "economy3")], 1, sum) - 3) * (10/15) # Re-scale to 0-10
D$economyIndex <- (D$economyIndex - weighted.mean(D$economyIndex[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$economyIndex[D$period == 0], D$weight.lang[D$period == 0]))

# Create cultural threat index / standardize
D$culture1 <- abs(D$culture1 - 7) # Reverse scale direction
D$cultureIndex <- (apply(D[, c("culture1", "culture2", "culture3")], 1, sum) - 3) * (10/15) # Re-scale to 0-10
D$cultureIndex <- abs(D$cultureIndex - 10) # Flip direction so higher values = greater threat
D$cultureIndex <- (D$cultureIndex - weighted.mean(D$cultureIndex[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$cultureIndex[D$period == 0], D$weight.lang[D$period == 0]))

# Create security threat index / standardize
D$securityIndex <- (apply(D[, c("security1", "security2", "security3")], 1, sum) - 3) * (10/15) # Re-scale to 0-10
D$securityIndex <- (D$securityIndex - weighted.mean(D$securityIndex[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$securityIndex[D$period == 0], D$weight.lang[D$period == 0]))

# Retrospective national economic evaluation / standardize
D$economyRetrospectiveCanada <- (D$economyRetrospectiveCanada - weighted.mean(D$economyRetrospectiveCanada[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$economyRetrospectiveCanada[D$period == 0], D$weight.lang[D$period == 0]))

# Retrospective personal economic evaluation / standardize
D$economyRetrospectiveSelf <- (D$economyRetrospectiveSelf - weighted.mean(D$economyRetrospectiveSelf[D$period == 0], D$weight.lang[D$period == 0], na.rm = TRUE)) / sqrt(wtd.var(D$economyRetrospectiveSelf[D$period == 0], D$weight.lang[D$period == 0]))

D$allowRefugeesFactor <- factor(D$allowRefugees, levels = 1:6)
D$allowRefugeesBinary <- as.numeric(D$allowRefugees >= 4)


# Factor analyses for anxiety index
factanal(na.omit(D[, c("anxiety", "upset", "worry", "fear")]), factors = 1)
# (as noted in main manuscript and in Online Appendix A and F)
cronbach(na.omit(D[, c("anxiety", "upset", "worry", "fear")]))

# Factor analyses for anxiety index including the anger indicator
factanal(na.omit(D[, c("anxiety", "upset", "worry", "fear", "anger")]), factors = 1)
# Cronbach's alpha for anxiety index incl. anger indicator
# (as noted in Online Appendix F, and fn. 4 of the Online Appendix)
cronbach(na.omit(D[, c("anxiety", "upset", "worry", "fear", "anger")]))

# Factor analyses for sympathy index
factanal(na.omit(D[, c("sympathy", "indifference", "compassion", "sadness")]), factors = 1)
# Cronbach's alpha for sympathy index
# (as noted in main manuscript and in Online Appendix A)
cronbach(na.omit(D[, c("sympathy", "indifference", "compassion", "sadness")]))

# Factor analyses for security index
factanal(na.omit(D[, c("security1", "security2", "security3")]), factors = 1)
# Cronbach's alpha for security index
# (as noted in main manuscript and in Online Appendix A)
cronbach(na.omit(D[, c("security1", "security2", "security3")]))

# Factor analyses for cultural threat index
factanal(na.omit(D[, c("culture1", "culture2", "culture3")]), factors = 1)
# Cronbach's alpha for cultural threat index
# (as noted in main manuscript and in Online Appendix A)
cronbach(na.omit(D[, c("culture1", "culture2", "culture3")]))

# Factor analyses for economic threat index
factanal(na.omit(D[, c("economy1", "economy2", "economy3")]), factors = 1)
# Cronbach's alpha for economic threat index
# (as noted in main manuscript and in Online Appendix A)
cronbach(na.omit(D[, c("economy1", "economy2", "economy3")]))


# Total sample size
nrow(D) # 18,634

# Number of respondents before the attacks
sum(D$period == 0) # 1,152

# Number of respondents immediately after the attacks
sum(D$period == 1) # 2,448

# Number of respondents after that
sum(D$period > 1) # 15,034

# Average number of respondents per day after the immediate post-attack period
sum(D$period > 1) / 18 # 18 days


# Sample weight for the "language of use" variable for the raw data.
# Respondents' probability of being invited to the survey differed
# between Wave 1 and subsequent waves (by design), based on respondents'
# preferred survey language (French & English).
# See Online Appendix A for details.
surveyDesign <- svydesign(id = ~ 1, weights = ~ weight.lang, data = D)



##########
########## BALANCE TESTS
##########

# Create balance check table for immediate pre- and post-attacks period
getCheck <- function(x, name) {
  x <- as.formula(paste0("~ ", x))
  mean1 <- svymean(x, design = subset(surveyDesign, period == 0), na.rm = TRUE)
  mean1 <- mean1[length(mean1)]
  var1 <- svyvar(x, design = subset(surveyDesign, period == 0), na.rm = TRUE)[1]

  mean2 <- svymean(x, design = subset(surveyDesign, period == 1), na.rm = TRUE)
  mean2 <- mean2[length(mean2)]
  var2 <- svyvar(x, design = subset(surveyDesign, period == 1), na.rm = TRUE)[1]

  normalizedDiff <- (mean1 - mean2) / sqrt((var1 + var2)/2)

  return(paste0(name, " & ", round(mean1, 2), " & ", round(mean2, 2), " & ",
                round(normalizedDiff, 2)))
}

# Table A1
cat(c(getCheck("selfPlacement", "Political ideology (0 = left, 10 = right)"),
      getCheck("gender == 'Women'", "Female"),
      getCheck("gender == 'Men'", "Male"),
      getCheck("ageGroup == '18-29'", "Age 18-29"),
      getCheck("ageGroup == '30-39'", "Age 30-39"),
      getCheck("ageGroup == '40-49'", "Age 40-49"),
      getCheck("ageGroup == '50-64'", "Age 50-64"),
      getCheck("ageGroup == '65+'", "Age 65+"),
      getCheck("education == 'High school or below'", "High school or below"),
      getCheck("education == 'College'", "College"),
      getCheck("education == 'University degree'", "University degree"),
      getCheck("region == 'Quebec'", "Quebec"),
      getCheck("region == 'West'", "West"),
      getCheck("region == 'Ontario'", "Ontario"),
      getCheck("region == 'Atlantic'", "Atlantic"),
      getCheck("motherTongue == 'English'", "English"),
      getCheck("motherTongue == 'French'", "French"),
      getCheck("motherTongue == 'Other language'", "Other language")
      ), sep = " \\\\\n")


# Balance checks for all periods (Table A3)
getFullCheck <- function(x, name) {
  x <- as.formula(paste0("~ ", x))
  mean1 <- sapply(0:10, function(p) {
      value = svymean(x, design = subset(surveyDesign, period == p), na.rm = TRUE)
      return(value[length(value)])
    })

  return(paste0(c(name, round(mean1, 2)), collapse = " & "))

}

# Table A3
cat(c(getFullCheck("selfPlacement", "Political ideology (0 = left, 10 = right)"),
      getFullCheck("gender == 'Women'", "Female"),
      getFullCheck("gender == 'Men'", "Male"),
      getFullCheck("ageGroup == '18-29'", "Age 18-29"),
      getFullCheck("ageGroup == '30-39'", "Age 30-39"),
      getFullCheck("ageGroup == '40-49'", "Age 40-49"),
      getFullCheck("ageGroup == '50-64'", "Age 50-64"),
      getFullCheck("ageGroup == '65+'", "Age 65+"),
      getFullCheck("education == 'High school or below'", "High school or below"),
      getFullCheck("education == 'College'", "College"),
      getFullCheck("education == 'University degree'", "University degree"),
      getFullCheck("region == 'Quebec'", "Quebec"),
      getFullCheck("region == 'West'", "West"),
      getFullCheck("region == 'Ontario'", "Ontario"),
      getFullCheck("region == 'Atlantic'", "Atlantic"),
      getFullCheck("motherTongue == 'English'", "English"),
      getFullCheck("motherTongue == 'French'", "French"),
      getFullCheck("motherTongue == 'Other language'", "Other language")
      ), sep = " \\\\\n")


# Is there differential non-response in the immediate pre- and post-Paris
# attacks period
model.response.1 <- glm(afterParis ~ gender +
                                     ageGroup +
                                     education +
                                     motherTongue +
                                     region +
                                     selfPlacement +
                                     languageOfUse,
                        data = subset(D, period %in% c(0, 1)),
                        family = binomial)
summary(model.response.1)

# Table A2
stargazer(model.response.1,
          title = "",
          style = "apsr",
          single.row = TRUE,
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c(""),
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Survey language: French",
                               "Intercept"))

# Likelihood ratio tests comparing pre-attack sample to post-attack samples
lr_test <- function(t) {
  model.response <- glm(afterParis ~ gender +
                                     ageGroup +
                                     education +
                                     motherTongue +
                                     region +
                                     selfPlacement +
                                     languageOfUse,
                        data = subset(D, period %in% c(0, t)),
                        family = binomial)

  model.response.null <- glm(afterParis ~ languageOfUse,
                             data = subset(D, period %in% c(0, t) &
                                           !is.na(gender) &
                                           !is.na(ageGroup) &
                                           !is.na(education) &
                                           !is.na(motherTongue) &
                                           !is.na(region) &
                                           !is.na(selfPlacement) &
                                           !is.na(languageOfUse)),
                             family = binomial)
  return(lrtest(model.response, model.response.null)$`Pr(>Chisq)`[2])
}

# As noted in the main article and in Online Appendix B
# Comparison of pre-attack to immediate post-attack period
lr_test(1)

# As noted in the main article and in Online Appendix B
# Comparison of pre-attack to complete post-attack period
lr_test(1:10)

# Online Appendix B
# Comparison of pre-attack to each post-attack period separately, correcting
# p-values for the fact of making 10 separate comparisons
p.adjust(c(lr_test(1), lr_test(2), lr_test(3), lr_test(4), lr_test(5),
           lr_test(6), lr_test(7), lr_test(8), lr_test(9), lr_test(10)),
         method = "bonferroni")

# Equivalence tests for balance
E0 <- subset(D, period == 0) # Pre-attack sample
E1 <- subset(D, period == 1) # Immediate post-attack sample

# The TOSTER library does not allow for the use of survey weights.
# In the immediate post-attack period, we therefore remove, at random,
# English-language speakers to create balance on this variable between the pre-
# and post-attacks sample. This was done by design in the survey out-go, and
# which we use weights for otherwise. This procedure lowers statistical power
# but allows for the use of unweighted testing in this section

# Number to sample to have equal proportions French-language and English-language
num_samp <- round(sum(E1$languageOfUse == "English") * mean(E0$languageOfUse == "English") / mean(E0$languageOfUse == "French") *
                  (sum(E1$languageOfUse == "French") / sum(E1$languageOfUse == "English")))

set.seed(0)
E1 <- rbind(E1[E1$languageOfUse == "French", ],
            E1[sample(which(E1$languageOfUse == "English"), num_samp), ])
E <- rbind(E0, E1)

# Use a standardized version of the ideology scale
E$selfPlacement_sd <- scale(E$selfPlacement)

G_E <- data.frame(variable = c("Political ideology (0 = left, 10 = right)",
                               "Female",
                               "Male",
                               "Age 18-29",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "High school or below",
                               "College",
                               "University degree",
                               "Quebec",
                               "West",
                               "Ontario",
                               "Atlantic",
                               "English",
                               "French",
                               "Other language"),
                  est = NA, lower = NA, upper = NA)

# Get the equivalence test and equivalence intervals for a specified variable
get_eqb <- function(data, name, var_name, level = NA, group = "period",
                    type = "prop", alpha = 0.05,
                    low_eqbound = -0.25, high_eqbound = 0.25) {

    est <- lower <- upper <- std_dev <- est_sd <- lower_sd <- upper_sd <- NA

    if(type == "t") {
        tost_est <- dataTOSTtwo(E, var_name, group, alpha = alpha,
                                low_eqbound = -0.25, high_eqbound = 0.25)
        std_dev <- sd(data[, var_name], na.rm = TRUE)
        est <- as.numeric(diff(t(tost_est$desc$asDF[, c("m[2]", "m[1]")])))
        lower <- as.numeric(tost_est$eqb$asDF["cil[raw]"])
        upper <- as.numeric(tost_est$eqb$asDF["ciu[raw]"])
        p_upper <- as.numeric(tost_est$tost$asDF["p[1]"])
        p_lower <- as.numeric(tost_est$tost$asDF["p[2]"])
    } else if(type == "prop") {
        tost_est <- datatosttwoprop(E, var_name, level, group, alpha, 
                                    low_eqbound = -0.25, high_eqbound = 0.25)
        std_dev <- sd(data[, var_name] == level, na.rm = TRUE)
        est <- as.numeric(diff(t(tost_est$desc$asDF[, c("prop[2]", "prop[1]")])))
        lower <- as.numeric(tost_est$eqb$asDF["cil"])
        upper <- as.numeric(tost_est$eqb$asDF["ciu"])
        p_upper <- as.numeric(tost_est$tost$asDF["p[1]"])
        p_lower <- as.numeric(tost_est$tost$asDF["p[2]"])
    }

    out <- data.frame(name = name, est = est, lower = lower, upper = upper,
                      std_dev = std_dev, est_sd = est / std_dev,
                      lower_sd = lower / std_dev, upper_sd = upper / std_dev,
                      p_upper = p_upper, p_lower = p_lower)
    row.names(out) <- 1:nrow(out)

    return(out)

}

# Get data for equivalence tests for graphing Figure A1
G_E <- rbind(get_eqb(E, "Political ideology (0 = left, 10 = right)", "selfPlacement", type = "t"),
             get_eqb(E, "Female", "gender", "Women"),
             get_eqb(E, "Male", "gender", "Men"),
             get_eqb(E, "Age 18-29", "ageGroup", "18-29"),
             get_eqb(E, "Age 30-39", "ageGroup", "30-39"),
             get_eqb(E, "Age 40-49", "ageGroup", "40-49"),
             get_eqb(E, "Age 50-64", "ageGroup", "50-64"),
             get_eqb(E, "Age 65+", "ageGroup", "65+"),
             get_eqb(E, "High school or below", "education", "High school or below"),
             get_eqb(E, "College", "education", "College"),
             get_eqb(E, "University degree", "education", "University degree"),
             get_eqb(E, "Quebec", "region", "Quebec"),
             get_eqb(E, "West", "region", "West"),
             get_eqb(E, "Ontario", "region", "Ontario"),
             get_eqb(E, "English", "motherTongue", "English"),
             get_eqb(E, "French", "motherTongue", "French"),
             get_eqb(E, "Other language", "motherTongue", "Other language"))

G_E$name <- factor(G_E$name, levels = rev(c("Political ideology (0 = left, 10 = right)",
                                        "Female", "Male",
                                        "Age 18-29", "Age 30-39", "Age 40-49", 
                                        "Age 50-64", "Age 65+",
                                        "High school or below",
                                        "College",
                                        "University degree",
                                        "Quebec", "West", "Ontario",
                                        "English", "French", "Other language")))

# Figure A1
pdf("Graphs/Figure_A1-Raw.pdf", 6.5, 3.75, useDingbats = FALSE)
ggplot(G_E, aes(x = est_sd, y = name)) +
  my.theme(borderless = 2, remove.ticks.y = TRUE,
           grid.y_colour = "grey80", grid.y_linetype = 3,
           base_size = 9) +
  coord_cartesian(xlim = c(-0.36, 0.36)) +
  scale_x_continuous(breaks = c(-0.36, -0.3, -0.25, -0.2, -0.1, 0,
                                0.1, 0.2, 0.25, 0.3, 0.36),
                     labels = c("-.36", "-.3", "-.25", "-.2", "-.1", "0",
                                ".1", ".2", ".25", ".3", ".36")) +
  labs(x = "Equivalence range (in std. dev.)", y = "") +
  geom_vline(xintercept = c(-0.36, -0.25, 0.25, 0.36), linetype = 3, size = 0.25) +
  geom_segment(aes(x = lower_sd, xend = upper_sd, yend = name),
               size = 4, color = "grey90") +
  geom_point() +
  geom_vline(xintercept = 0, linetype = 5, size = 0.25)
dev.off()

# p values testing whether imbalance on each variable is significantly less than
# epsilon = 0.25 SD, and greater than epsilon = 0.25 SD
# (as noted in Online Appendix B)
all(G_E[, c("p_lower", "p_upper")] < 0.001)


##########
########## GENERATE TIME SERIES DATA FOR GRAPHING
##########

# Set up the time series data for graphing
G <- data.frame(day = c(-1, 1, 3, 5, 7, 9, 11, 13, 15, 17, 19),
                treatment = c(rep("Before attack", 1), rep("After attack", 10)),
                sampleSize = c(length(na.omit(D$allowRefugees[D$period == 0 & D$afterParis == 0])),
                               length(na.omit(D$allowRefugees[D$period == 1 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 2 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 3 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 4 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 5 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 6 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 7 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 8 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 9 & D$afterParis == 1])),
                               length(na.omit(D$allowRefugees[D$period == 10 & D$afterParis == 1]))),
                allowRefugees = c(summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                  summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                allowRefugeesSE = c(summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                    summary(svyglm(allowRefugees ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                contactMP = c(summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                              summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                contactMPSE = c(summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                summary(svyglm(contactMP ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                donateIntegrate = c(summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                    summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                donateIntegrateSE = c(summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                      summary(svyglm(donateIntegrate ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                anxietyIndex = c(summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                anxietyIndexSE = c(summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                anxietyWithAnger = c(summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                 summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                anxietyWithAngerSE = c(summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                   summary(svyglm(anxietyWithAnger ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                angerStd = c(summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                 summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                angerStdSE = c(summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                   summary(svyglm(angerStd ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                sympathyIndex = c(summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                  summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                sympathyIndexSE = c(summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                    summary(svyglm(sympathyIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                securityIndex = c(summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                  summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                securityIndexSE = c(summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                    summary(svyglm(securityIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                economyIndex = c(summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                 summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                economyIndexSE = c(summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                   summary(svyglm(economyIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                cultureIndex = c(summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                 summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                cultureIndexSE = c(summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                   summary(svyglm(cultureIndex ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                economyRetrospectiveSelf = c(summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                             summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                economyRetrospectiveSelfSE = c(summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                               summary(svyglm(economyRetrospectiveSelf ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]),
                economyRetrospectiveCanada = c(summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[1],
                                               summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[1]),
                economyRetrospectiveCanadaSE = c(summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 0 & afterParis == 0)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 1 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 2 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 3 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 4 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 5 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 6 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 7 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 8 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 9 & afterParis == 1)))$coef[2],
                                                 summary(svyglm(economyRetrospectiveCanada ~ 1, design = subset(surveyDesign, period == 10 & afterParis == 1)))$coef[2]))

G <- melt(G, id.vars = c("day", "sampleSize", "treatment",
                         "allowRefugeesSE",
                         "contactMPSE",
                         "donateIntegrateSE",
                         "anxietyIndexSE",
                         "anxietyWithAngerSE",
                         "angerStdSE",
                         "sympathyIndexSE",
                         "securityIndexSE",
                         "economyIndexSE",
                         "cultureIndexSE",
                         "economyRetrospectiveCanadaSE",
                         "economyRetrospectiveSelfSE"))

# Rename variables to appear as desired in figure facets/headers
G$variable <- as.character(G$variable)
G$variable[G$variable == "allowRefugees"] <- "Support for refugee resettlement"
G$variable[G$variable == "contactMP"] <- "Willingness to contact MP\nregarding refugee resettlement"
G$variable[G$variable == "donateIntegrate"] <- "Would consider donating to help resettlement"
G$variable[G$variable == "anxietyIndex"] <- "Anxiety"
G$variable[G$variable == "anxietyWithAnger"] <- "Anxiety (incl. anger)"
G$variable[G$variable == "angerStd"] <- "Anger"
G$variable[G$variable == "sympathyIndex"] <- "Sympathy"
G$variable[G$variable == "securityIndex"] <- "Security Threat"
G$variable[G$variable == "economyIndex"] <- "Economic Threat"
G$variable[G$variable == "cultureIndex"] <- "Cultural Threat"
G$variable[G$variable == "economyRetrospectiveSelf"] <- "Retrospective evaluation\nof personal financial situation"
G$variable[G$variable == "economyRetrospectiveCanada"] <- "Retrospective evaluation\nof national economy"

# Change the factor levels to order as desired in the faceted figures
G$variable <- factor(G$variable,
                     levels = c("Support for refugee resettlement",
                                "Willingness to contact MP\nregarding refugee resettlement",
                                "Would consider donating to help resettlement",
                                "Anxiety",
                                "Anxiety (incl. anger)",
                                "Anger",
                                "Sympathy",
                                "Security Threat",
                                "Cultural Threat",
                                "Economic Threat",
                                "Retrospective evaluation\nof personal financial situation",
                                "Retrospective evaluation\nof national economy"))


##########
########## MAIN REGRESSION MODELS
##########

##### ANXIETY
# Regression Table A4 Model (1)
model.anxiety <- glm(anxietyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.anxiety)

# Regression Table A7 Model (2)
model.anxiety.immediate.responders <- glm(anxietyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
# summary(model.anxiety.immediate.responders)

# Regression Table A14 Model (3) - Factorial experiment (interaction)
model.anxiety.factorial1 <- glm(anxietyIndex ~ afterParis * treatment1 + afterParis * treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.anxiety.factorial1)

# Regression Table A14 Model (2) - Factorial experiment (additive)
model.anxiety.factorial2 <- glm(anxietyIndex ~ afterParis + treatment1 + treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.anxiety.factorial2)

# Regression Table A10 Model (2)
model.anxiety.anger <- glm(anxietyWithAnger ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.anxiety.anger)

# Regression Table A10 Model (3)
model.angerStd <- glm(angerStd ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.angerStd)


##### SYMPATHY
# Regression Table A4 Model (2)
model.sympathy <- glm(sympathyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.sympathy)

# Regression Table A7 Model (4)
model.sympathy.immediate.responders <- glm(sympathyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
# summary(model.sympathy.immediate.responders)

# Regression Table A15 Model (3) - Factorial experiment (interaction)
model.sympathy.factorial1 <- glm(sympathyIndex ~ afterParis*treatment1 + afterParis*treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.sympathy.factorial1)

# Regression Table A15 Model (2) - Factorial experiment (additive)
model.sympathy.factorial2 <- glm(sympathyIndex ~ afterParis + treatment1 + treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.sympathy.factorial2)


##### SECURITY THREAT
# Regression Table A4 Model (3)
model.security <- glm(securityIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.security)

# Regression Table A7 Model (6)
model.security.immediate.responders <- glm(securityIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
# summary(model.security.immediate.responders)

# Regression Table A16 Model (3) - Factorial experiment (interaction)
model.security.factorial1 <- glm(securityIndex ~ afterParis*treatment1 + afterParis*treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.security.factorial1)

# Regression Table A16 Model (3) - Factorial experiment (additive)
model.security.factorial2 <- glm(securityIndex ~ afterParis + treatment1 + treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.security.factorial2)


##### CULTURAL THREAT
# Regression Table A4 Model (4)
model.culture <- glm(cultureIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.culture)

# Regression Table A7 Model (8)
model.culture.immediate.responders <- glm(cultureIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
# summary(model.culture.immediate.responders)

# Regression Table A17 Model (3) - Factorial experiment (interaction)
model.culture.factorial1 <- glm(cultureIndex ~ afterParis*treatment1 + afterParis*treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.culture.factorial1)

# Regression Table A17 Model (3) - Factorial experiment (interaction)
model.culture.factorial2 <- glm(cultureIndex ~ afterParis + treatment1 + treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.culture.factorial2)


##### ECONOMIC THREAT
# Regression Table A13 Model (5)
model.economy <- glm(economyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.economy)

# Regression Table A7 Model (10)
model.economy.immediate.responders <- glm(economyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
# summary(model.economy.immediate.responders)

# Not shown
model.economy.factorial1 <- glm(economyIndex ~ afterParis*treatment1 + afterParis*treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.economy.factorial1)

# Not shown
model.economy.factorial2 <- glm(economyIndex ~ afterParis + treatment1 + treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement, data = subset(D, period %in% c(0, 1)))
# summary(model.economy.factorial2)


# Table A4
stargazer(model.anxiety, model.sympathy, model.security, model.culture,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = TRUE,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Intercept"))
# coefficient and p-value for anxiety as noted in main article
round(coef(summary(model.anxiety))["afterParis", c(1, 4)], 4)
# coefficient and p-value for sympathy as noted in main article
round(coef(summary(model.sympathy))["afterParis", c(1, 4)], 4)
# coefficient and p-value for security as noted in main article
round(coef(summary(model.security))["afterParis", c(1, 4)], 4)
# coefficient and p-value for culture as noted in main article
round(coef(summary(model.culture))["afterParis", c(1, 4)], 4)



# Table A13 (includes economic threat)
stargazer(model.anxiety, model.sympathy, model.security,
          model.culture, model.economy,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = TRUE,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Intercept"))


# Table A7
stargazer(model.anxiety,
          model.anxiety.immediate.responders,
          model.sympathy,
          model.sympathy.immediate.responders,
          model.security,
          model.security.immediate.responders,
          model.culture,
          model.culture.immediate.responders,
          model.economy,
          model.economy.immediate.responders,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = TRUE,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Intercept"))


# Table A10
stargazer(model.anxiety, model.anxiety.anger, model.angerStd,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = TRUE,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Intercept"))


# Set up the max and min y-axis limits for Figure 2
G$max <- G$min <- as.numeric(NA)
G$min[G$variable == "Anxiety"] <- -0.05
G$max[G$variable == "Anxiety"] <- 0.35

G$min[G$variable == "Anxiety (incl. anger)"] <- -0.05
G$max[G$variable == "Anxiety (incl. anger)"] <- 0.35

G$min[G$variable == "Sympathy"] <- -0.3
G$max[G$variable == "Sympathy"] <- 0.3

G$min[G$variable == "Security Threat"] <- -0.19
G$max[G$variable == "Security Threat"] <- 0.46

G$min[G$variable == "Cultural Threat"] <- -0.19
G$max[G$variable == "Cultural Threat"] <- 0.46

G$min[G$variable == "Economic Threat"] <- -0.35
G$max[G$variable == "Economic Threat"] <- 0.35

# Note: Most figures below are mildly edited in Adobe Illustrator after
#       being written out to pdf

# Figure 1
pdf("Graphs/Figure_1-Raw.pdf", 4.6, 2.1, useDingbats = FALSE)
ggplot(subset(G, variable %in% c("Anxiety", "Sympathy")), aes(x = day, y = value, size = treatment, colour = treatment, fill = treatment)) +
  my.theme(bordersize = 1, base_size = 9) +
  facet_wrap(~ variable, scales = "free", ncol = 2) +
  labs(x = "Days from attack", y = "Index value (std. units)", fill = "", colour = "") +
  # coord_cartesian(xlim = c(-2.55, 20.10), expand = FALSE) +
  scale_x_continuous(limits = c(-2.55, 20.10), expand = c(0, 0),
                     breaks = c(0, 5, 10, 15, 20, 25, 30),
                     labels = c("Paris\nattacks", "5", "10", "15", "20", "25", "30")) +
  scale_y_continuous(breaks = seq(-0.4, 0.4, by = 0.1)) +
  geom_hline(yintercept = 0, size = 0.25, linetype = 3, color = "grey80") +
  geom_vline(xintercept = 0, size = 0.25, linetype = 2, colour = "black") +
  geom_smooth(size = 3.5, fullrange = FALSE, colour = "grey90", method = "lm", formula = y ~ x + I(x^2), se = FALSE) +
  geom_linerange(data = subset(G, variable == "Anxiety"), aes(ymax = value + anxietyIndexSE * qnorm(0.95), ymin = value - anxietyIndexSE * qnorm(0.95)), size = 0.3) +
  geom_linerange(data = subset(G, variable == "Sympathy"), aes(ymax = value + sympathyIndexSE * qnorm(0.95), ymin = value - sympathyIndexSE * qnorm(0.95)), size = 0.3) +
  geom_point(shape = 21, size = 2.5) +
  geom_blank(aes(y = max)) +
  geom_blank(aes(y = min)) +
  scale_size_manual(values = c("Before attack" = 2.5, "After attack" = 2.25)) +
  scale_fill_manual(values = c("Before attack" = "white", "After attack" = "black")) +
  scale_colour_manual(values = c("Before attack" = "black", "After attack" = "black")) +
  theme(legend.position = "none",
        panel.spacing = unit(0.4, "lines"),
        plot.margin = unit(c(0, 0.35, 0.25, 0.25), "lines"))
dev.off()


# Figure 2
pdf("Graphs/Figure_2-Raw.pdf", 4.6, 2.1, useDingbats = FALSE)
ggplot(subset(G, variable %in% c("Security Threat", "Cultural Threat")), aes(x = day, y = value, size = treatment, colour = treatment, fill = treatment)) +
  my.theme(bordersize = 1, base_size = 9) +
  facet_wrap(~ variable, scales = "free_y", ncol = 3) +
  labs(x = "Days from attack", y = "Index value (std. units)", fill = "", colour = "") +
  scale_x_continuous(limits = c(-2.55, 20.10), expand = c(0, 0),
                     breaks = c(0, 5, 10, 15, 20, 25, 30),
                     labels = c("Paris\nattacks", "5", "10", "15", "20", "25", "30")) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(-0.5, 0.5, by = 0.1)) +
  geom_hline(yintercept = 0, size = 0.25, linetype = 3, color = "grey80") +
  geom_vline(xintercept = 0, size = 0.25, linetype = 2, colour = "black") +
  geom_smooth(size = 2.75, fullrange = FALSE, colour = "grey90", method = "lm", formula = y ~ x + I(x^2), se = FALSE) +
  geom_linerange(data = subset(G, variable == "Security Threat"), aes(ymax = value + securityIndexSE * qnorm(0.95), ymin = value - securityIndexSE * qnorm(0.95)), size = 0.3) +
  geom_linerange(data = subset(G, variable == "Cultural Threat"), aes(ymax = value + cultureIndexSE * qnorm(0.95), ymin = value - cultureIndexSE * qnorm(0.95)), size = 0.3) +
  # geom_linerange(data = subset(G, variable == "Economic Threat"), aes(ymax = value + economyIndexSE * qnorm(0.95), ymin = value - economyIndexSE * qnorm(0.95)), size = 0.3) +
  geom_point(shape = 21, size = 2.5) +
  geom_blank(aes(y = max)) +
  geom_blank(aes(y = min)) +
  scale_size_manual(values = c("Before attack" = 2.5, "After attack" = 2.25)) +
  scale_fill_manual(values = c("Before attack" = "white", "After attack" = "black")) +
  scale_colour_manual(values = c("Before attack" = "black", "After attack" = "black")) +
  theme(legend.position = "none",
        panel.spacing = unit(0.4, "lines"),
        plot.margin = unit(c(0, 0.35, 0.25, 0.25), "lines"))
dev.off()


# Figure A4
pdf("Graphs/Figure_A4-Raw.pdf", 3.5, 2.5, useDingbats = FALSE)
ggplot(subset(G, variable %in% c("Economic Threat")), aes(x = day, y = value, size = treatment, colour = treatment, fill = treatment)) +
  my.theme(bordersize = 1, base_size = 10) +
  facet_wrap(~ variable, scales = "free_y", ncol = 3) +
  labs(x = "Days from attack", y = "Index value (std. units)", fill = "", colour = "") +
  scale_x_continuous(limits = c(-2.55, 20.10), expand = c(0, 0), breaks = c(0, 5, 10, 15, 20, 25, 30)) +
  scale_y_continuous(breaks = seq(-0.5, 0.5, by = 0.1)) +
  geom_vline(xintercept = 0, size = 0.25, linetype = 2, colour = "black") +
  geom_smooth(size = 2.75, fullrange = FALSE, colour = "grey90", method = "lm", formula = y ~ x + I(x^2), se = FALSE) +
  # geom_linerange(data = subset(G, variable == "Security Threat"), aes(ymax = value + securityIndexSE * qnorm(0.95), ymin = value - securityIndexSE * qnorm(0.95)), size = 0.3) +
  # geom_linerange(data = subset(G, variable == "Cultural Threat"), aes(ymax = value + cultureIndexSE * qnorm(0.95), ymin = value - cultureIndexSE * qnorm(0.95)), size = 0.3) +
  geom_linerange(data = subset(G, variable == "Economic Threat"), aes(ymax = value + economyIndexSE * qnorm(0.95), ymin = value - economyIndexSE * qnorm(0.95)), size = 0.3) +
  geom_point(shape = 21, size = 2.5) +
  geom_blank(aes(y = max)) +
  geom_blank(aes(y = min)) +
  scale_size_manual(values = c("Before attack" = 2.5, "After attack" = 2.25)) +
  scale_fill_manual(values = c("Before attack" = "white", "After attack" = "black")) +
  scale_colour_manual(values = c("Before attack" = "black", "After attack" = "black")) +
  theme(legend.position = "none",
        plot.margin = unit(c(0.5, 0.35, 0.25, 0.25), "lines"))
dev.off()


# Figure 3 Panel A
# Note: Figures 3 is post-processed in Illustrator to compile
# panels A and B side-by-side, in addition to further labelling
pdf("Graphs/Figure_3A.pdf", 3.25, 2.4, useDingbats = FALSE)
ggplot(subset(G, variable %in% c("Support for refugee resettlement")), aes(x = day, y = value, size = treatment, colour = treatment, fill = treatment)) +
  my.theme(bordersize = 1, base_size = 10) +
  # facet_wrap(~ variable, scales = "free", ncol = 1) +
  labs(x = "Days from attack", y = "Support for refugees", fill = "", colour = "") +
  coord_cartesian(xlim = c(-2.55, 20.10), ylim = c(4.64, 5.21), expand = FALSE) +
  scale_x_continuous(breaks = c(0, 5, 10, 15, 20, 25, 30), labels = c("Paris\nattacks", "5", "10", "15", "20", "25", "30")) +
  scale_y_continuous(breaks = seq(4.7, 6, by = 0.1)) +
  geom_hline(yintercept = 0, size = 0.25, linetype = 3, color = "grey80") +
  geom_vline(xintercept = 0, size = 0.25, linetype = 2, colour = "black") +
  geom_smooth(size = 3, fullrange = FALSE, colour = "grey90", method = "lm", formula = y ~ x + I(x^2), se = FALSE) +
  geom_linerange(aes(ymax = value + allowRefugeesSE * qnorm(0.95), ymin = value - allowRefugeesSE * qnorm(0.95)), size = 0.3) +
  geom_point(shape = 21, size = 2.5) +
  scale_size_manual(values = c("Before attack" = 2.5, "After attack" = 2.25)) +
  scale_fill_manual(values = c("Before attack" = "white", "After attack" = "black")) +
  scale_colour_manual(values = c("Before attack" = "black", "After attack" = "black")) +
  theme(legend.position = "none",
        plot.margin = unit(c(0.5, 0.35, 0.25, 0.25), "lines"))
dev.off()


# Table A5 Model (1)
model.refugees1 <- polr(as.factor(allowRefugees) ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                        data = subset(D, period %in% c(0, 1)))
summary(model.refugees1, digits = 2)

# Table A8 Model (2)
model.refugees1.immediate.respondents <- polr(as.factor(allowRefugees) ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                                              data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
summary(model.refugees1.immediate.respondents, digits = 2)

# Table A18 Model (3) (interaction)
model.refugees1.factorial1 <- polr(as.factor(allowRefugees) ~ afterParis*treatment1 + afterParis*treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement,
                                   data = subset(D, period %in% c(0, 1)))
summary(model.refugees1.factorial1, digits = 2)

# Table A18 Model (2) (additive)
model.refugees1.factorial2 <- polr(as.factor(allowRefugees) ~ afterParis + treatment1 + treatment2 + gender + ageGroup + education + motherTongue + region + selfPlacement,
                                   data = subset(D, period %in% c(0, 1)))
summary(model.refugees1.factorial2, digits = 2)


# Table A5
stargazer(model.refugees1,
          single.row = TRUE,
          title = "Support for refugee resettlement before and after Paris terrorist attacks",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = FALSE,
          model.names = TRUE,
          ord.intercepts = TRUE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "\\midrule$\\kappa_1$",
                               "$\\kappa_2$",
                               "$\\kappa_3$",
                               "$\\kappa_4$",
                               "$\\kappa_5$"))


# Table A8
stargazer(model.refugees1, model.refugees1.immediate.respondents,
          # single.row = TRUE,
          title = "Support for refugee resettlement before and after Paris terrorist attacks",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = FALSE,
          model.names = TRUE,
          ord.intercepts = TRUE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "\\midrule$\\kappa_1$",
                               "$\\kappa_2$",
                               "$\\kappa_3$",
                               "$\\kappa_4$",
                               "$\\kappa_5$"))


# Point estimates for 6-category support for refugee resettlement outcome
G2 <- data.frame(category1 = c("Strongly\ndisagree",
                               "Somewhat\ndisagree",
                               "Slightly\ndisagree",
                               "Slightly\nagree",
                               "Somewhat\nagree",
                               "Strongly\nagree"),
                 category2 = rep(c("Disagree", "Agree"), each = 3),
                 value = NA,
                 lower = NA,
                 upper = NA)

get_sim <- function(fitted_model, n_sims) {

  # Simulate coefficients
  set.seed(1)
  coefs_sim <- rmvnorm(n_sims, c(fitted_model$coefficients, fitted_model$zeta), vcov(fitted_model))

  # Cut points
  zeta_sim <- coefs_sim[, (ncol(coefs_sim) - length(fitted_model$zeta) + 1):ncol(coefs_sim)]

  # Linear predictor
  coefs_sim <- coefs_sim[, 1:(ncol(coefs_sim) - length(fitted_model$zeta))]

  # Model matrices with predictors for before and after the attacks
  data_sim_0 <- data_sim_1 <- subset(D, !is.na(weight))
  data_sim_0$afterParis <- 0
  data_sim_1$afterParis <- 1
  weight_sim <- data_sim_0$weight

  sims <- t(sapply(1:n_sims, function(x) {
              fitted_model$coefficients <- coefs_sim[x, ]
              fitted_model$zeta <- zeta_sim[x, ]
              apply(predict(fitted_model,
                               newdata = data_sim_1,
                               type = "prob",
                               na.action = "na.pass") -
                    predict(fitted_model,
                            newdata = data_sim_0,
                            type = "prob",
                            na.action = "na.pass"), 2,
                    function(x) weighted.mean(x, weight_sim, na.rm = TRUE)) }))

  return(data.frame(value = colMeans(sims),
                    lower = apply(sims, 2, function(x) quantile(x, 0.05)),
                    upper = apply(sims, 2, function(x) quantile(x, 0.95))))

}

# Get predicted probabilities for the 6-item refugee resettlement outcome
# for before and after the attacks, simulating parameters from the model
# Roughly 8 minutes to run
G2[, c("value", "lower", "upper")] <- get_sim(model.refugees1, 1000) * 100

G2$category1 <- factor(G2$category1,
                       levels = c("Strongly\ndisagree",
                                  "Somewhat\ndisagree",
                                  "Slightly\ndisagree",
                                  "Slightly\nagree",
                                  "Somewhat\nagree",
                                  "Strongly\nagree"))

# Magnitude of the effect of the attacks on support for refugee resettlement
# as noted in the text of the main article, and in Panel B of Figure 3
sum(G2$value[4:6])
# p-value as noted in the main article
pt(coef(summary(model.refugees1))["afterParis", 3],
   df = df.residual(model.refugees1))


# Figure 3 Panel B
# Note: Figures 3 is post-processed in Illustrator to compile
# panels A and B side-by-side, in addition to further labelling
pdf("Graphs/Figure_3B.pdf", 3.25, 2.4, useDingbats = FALSE)
ggplot(G2, aes(x = category1, y = value, fill = category2)) +
  my.theme(borderless = 2, base_size = 9) +
  labs(color = "", fill = "", y = "Effect size (%-points)", x = "") +
  coord_cartesian(ylim = c(-8, 6), expand = FALSE) +
  scale_y_continuous(breaks = c()) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.7, colour = "black", size = 0.25) +
  geom_linerange(aes(ymax = upper, ymin = lower), size = 0.5,
               position = position_dodge(width = 0.9), color = "black") +
  geom_hline(yintercept = 0, size = 0.25) +
  geom_text(aes(y = upper + 0.75, label = round(value)), position = position_dodge(width = 0.9), colour = "black", size = 3.2) +
  scale_fill_manual(values = c("Disagree" = "grey70",
                               "Agree" = "white")) +
  theme(legend.position = "none",
        legend.direction = "vertical",
        plot.margin = unit(c(-0.75, 0.25, -0.25, 0.25), "lines"))
dev.off()




##########
########## CAUSAL MECHANISMS ANALYSIS
##########

set.seed(1)
n_sims <- 5000

##### Anxiety causal mediation regressions

# Mediator regression
model.anxiety <- lm(anxietyIndex ~ afterParis +
                                   gender + ageGroup + education + motherTongue + region + selfPlacement,
                    data = subset(D, allSocio == 1 & !is.na(anxietyIndex) & !is.na(sympathyIndex) &
                                     !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
# Outcome regression
model.anxiety.outcome <- polr(allowRefugeesFactor ~ afterParis + anxietyIndex +
                                                    gender + ageGroup + education + motherTongue + region + selfPlacement,
                              data = subset(D, allSocio == 1 & !is.na(anxietyIndex) & !is.na(sympathyIndex) &
                                               !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.anxiety.outcome)

# Mediation (30 minutes to run; load pre-run version below)
mediate.anxiety <- mediate(model.anxiety, model.anxiety.outcome, boot = TRUE,
                           treat = "afterParis", mediator = "anxietyIndex",
                           sims = n_sims)
# saveRDS(mediate.anxiety, "Models/mediate.anxiety.rds")
# mediate.anxiety <- readRDS("Models/mediate.anxiety.rds")

# Simulation results (4:6 represent Slightly/Somewhat/Strongly agree)
# We thus averaging over these for z = 0 and z = 1
sims.anxiety.acme <- rbind(mediate.anxiety$d0.sims[, 4:6],
                           mediate.anxiety$d1.sims[, 4:6])
sims.anxiety.ade <- rbind(mediate.anxiety$z0.sims[, 4:6],
                          mediate.anxiety$z1.sims[, 4:6])
sims.anxiety.tau <- rbind(mediate.anxiety$tau.sims[, 4:6])

# Table A12 (Anxiety ACME)
# As noted in Online Appendix G
mean(apply(sims.anxiety.acme, 1, sum)) # ACME
quantile(apply(sims.anxiety.acme, 1, sum), c(0.025, 0.975)) # ACME CIs

# ADE (Not shown)
mean(apply(sims.anxiety.ade, 1, sum))
quantile(apply(sims.anxiety.ade, 1, sum), c(0.025, 0.975))

# Total effect (Not shown)
mean(apply(sims.anxiety.tau, 1, sum))
quantile(apply(sims.anxiety.tau, 1, sum), c(0.025, 0.975))


# Sensitivity analysis
model.anxiety.outcome.ols <- lm(allowRefugees ~ afterParis + anxietyIndex +
                                                gender + ageGroup + education + motherTongue + region + selfPlacement,
                                data = subset(D, allSocio == 1 & !is.na(anxietyIndex) & !is.na(sympathyIndex) &
                                               !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.anxiety.outcome.ols)

# Mediation
mediate.anxiety.ols <- mediate(model.anxiety, model.anxiety.outcome.ols, boot = FALSE, treat = "afterParis",
                               mediator = "anxietyIndex", sims = n_sims)

mediate.anxiety.ols.sens <- medsens(mediate.anxiety.ols, rho.by = 0.01,
                                    effect.type = "indirect")
# Correlation (with unobserved confounder) at which the AMCE is zero
# (as noted in Online Appendix G)
summary(mediate.anxiety.ols.sens) # "Rho at which ACME = 0: -0.62"


# Figure A3 Panel Top Left
# Note: Panels in A3 are assembled in Illustrator
pdf("Graphs/Figure_A3_Top_Left-Raw.pdf", 3.5, 3.9, useDingbats = FALSE) # Shrink to 3.25 height,
plot(mediate.anxiety.ols.sens, main = "ACME: Anxiety",
     cex.main = 0.9, cex.axis = 0.85)
dev.off()


##### Sympathy causal mediation regressions

# Mediator regression
model.sympathy <- lm(sympathyIndex ~ afterParis +
                                     gender + ageGroup + education + motherTongue + region + selfPlacement,
                     data = subset(D, allSocio == 1 & !is.na(anxietyIndex) & !is.na(sympathyIndex) &
                                      !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
# Outcome regression
model.sympathy.outcome <- polr(allowRefugeesFactor ~ afterParis + sympathyIndex +
                                                     gender + ageGroup + education + motherTongue + region + selfPlacement,
                               data = subset(D, allSocio == 1 & !is.na(anxietyIndex) & !is.na(sympathyIndex) &
                                                !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.sympathy.outcome)

# Mediation (30 minutes to run; load pre-run version below)
mediate.sympathy <- mediate(model.sympathy, model.sympathy.outcome, boot = TRUE, treat = "afterParis",
                            mediator = "sympathyIndex", sims = n_sims)
# saveRDS(mediate.sympathy, "Models/mediate.sympathy.rds")
# mediate.sympathy <- readRDS("Models/mediate.sympathy.rds")

# Simulation results (4:6 represent Slightly/Somewhat/Strongly agree)
# Averaging over z = 0 and z =1
sims.sympathy.acme <- rbind(mediate.sympathy$d0.sims[, 4:6],
                            mediate.sympathy$d1.sims[, 4:6])
sims.sympathy.ade <- rbind(mediate.sympathy$z0.sims[, 4:6],
                           mediate.sympathy$z1.sims[, 4:6])
sims.sympathy.tau <- rbind(mediate.sympathy$tau.sims[, 4:6])


# Table A12 (Sympathy ACME)
# As noted in Online Appendix G
mean(apply(sims.sympathy.acme, 1, sum)) # ACME
quantile(apply(sims.sympathy.acme, 1, sum), c(0.025, 0.975)) # ACME CIs

# ADE
mean(apply(sims.sympathy.ade, 1, sum))
quantile(apply(sims.sympathy.ade, 1, sum), c(0.025, 0.975))

# Total effect
mean(apply(sims.sympathy.tau, 1, sum))
quantile(apply(sims.sympathy.tau, 1, sum), c(0.025, 0.975))


# Sensitivity analysis
model.sympathy.outcome.ols <- lm(allowRefugees ~ afterParis + sympathyIndex +
                                                 gender + ageGroup + education + motherTongue + region + selfPlacement,
                                 data = subset(D, allSocio == 1 & !is.na(anxietyIndex) & !is.na(sympathyIndex) &
                                               !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.sympathy.outcome.ols)
# Mediation
mediate.sympathy.ols <- mediate(model.sympathy, model.sympathy.outcome.ols, boot = FALSE, treat = "afterParis",
                                mediator = "sympathyIndex", sims = n_sims)

mediate.sympathy.ols.sens <- medsens(mediate.sympathy.ols, rho.by = 0.01,
                                     effect.type = "indirect")
summary(mediate.sympathy.ols.sens)

# Figure A3 Panel Top Right
# Note: Panels in A3 are assembled in Illustrator
pdf("Graphs/Figure_A3_Top_Right-Raw.pdf", 3.5, 3.9, useDingbats = FALSE) # Shrink to 3.25 height
plot(mediate.sympathy.ols.sens, main = "Sympathy",
     cex.main = 0.9, cex.axis = 0.85)
dev.off()


##### Security threat causal mediation regressions

# Mediator regression
model.security <- lm(securityIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                     data = subset(D, allSocio == 1 & !is.na(allowRefugeesFactor) & period %in% c(0, 1)))

# Outcome regression
model.security.outcome <- polr(allowRefugeesFactor ~ afterParis + securityIndex +
                                                     gender + ageGroup + education + motherTongue + region + selfPlacement,
                               data = subset(D, allSocio == 1 & !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.security.outcome)

# Mediation (30 minutes to run; load pre-run version below)
mediate.security <- mediate(model.security, model.security.outcome, boot = TRUE, treat = "afterParis",
                            mediator = "securityIndex", sims = n_sims)
# saveRDS(mediate.security, "Models/mediate.security.rds")
# mediate.security <- readRDS("Models/mediate.security.rds")


# Simulation results (4:6 represent Slightly/Somewhat/Strongly agree)
# Averaging over z = 0 and z =1
sims.security.acme <- rbind(mediate.security$d0.sims[, 4:6],
                            mediate.security$d1.sims[, 4:6])
sims.security.ade <- rbind(mediate.security$z0.sims[, 4:6],
                           mediate.security$z1.sims[, 4:6])
sims.security.tau <- rbind(mediate.security$tau.sims[, 4:6])


# Table A12 (Security threat ACME)
# As noted in Online Appendix G
mean(apply(sims.security.acme, 1, sum)) # ACME
quantile(apply(sims.security.acme, 1, sum), c(0.025, 0.975)) # ACME CIs

# ADE
mean(apply(sims.security.ade, 1, sum))
quantile(apply(sims.security.ade, 1, sum), c(0.025, 0.975))

# Total effect
mean(apply(sims.security.tau, 1, sum))
quantile(apply(sims.security.tau, 1, sum), c(0.025, 0.975))


# Sensitivity analysis
model.security.outcome.ols <- lm(allowRefugees ~ afterParis + securityIndex +
                                                 gender + ageGroup + education + motherTongue + region + selfPlacement,
                                 data = subset(D, allSocio == 1 & !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.security.outcome.ols)

# Mediation
mediate.security.ols <- mediate(model.security, model.security.outcome.ols, boot = FALSE, treat = "afterParis",
                                mediator = "securityIndex", sims = n_sims)

mediate.security.ols.sens <- medsens(mediate.security.ols, rho.by = 0.01,
                                     effect.type = "indirect")
# Correlation (with unobserved confounder) at which the AMCE is zero
# (as noted in Online Appendix G)
summary(mediate.security.ols.sens) # "Rho at which ACME = 0: -0.71"

# Figure A3 Panel Bottom Left
# Note: Panels in A3 are assembled in Illustrator
pdf("Graphs/Figure_A3_Bottom_Left-Raw.pdf", 3.5, 3.9, useDingbats = FALSE) # Shrink to 3.25 height
plot(mediate.security.ols.sens, main = "Security",
     cex.main = 0.9, cex.axis = 0.85)
dev.off()


##### Cultural threat causal mediation regressions

# Mediator regression
model.culture <- lm(cultureIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                    data = subset(D, allSocio == 1 & !is.na(allowRefugeesFactor) & period %in% c(0, 1)))

# Outcome regression
model.culture.outcome <- polr(allowRefugeesFactor ~ afterParis + cultureIndex +
                                                    gender + ageGroup + education + motherTongue + region + selfPlacement,
                              data = subset(D, allSocio == 1 & !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.culture.outcome)

# Mediation (30 minutes to run; load pre-run version below)
mediate.culture <- mediate(model.culture, model.culture.outcome, boot = TRUE, treat = "afterParis",
                           mediator = "cultureIndex", sims = n_sims)
# saveRDS(mediate.culture, "Models/mediate.culture.rds")
# mediate.culture <- readRDS("Models/mediate.culture.rds")


# Simulation results (4:6 represent Slightly/Somewhat/Strongly agree)
# Averaging over z = 0 and z =1
sims.culture.acme <- rbind(mediate.culture$d0.sims[, 4:6],
                           mediate.culture$d1.sims[, 4:6])
sims.culture.ade <- rbind(mediate.culture$z0.sims[, 4:6],
                          mediate.culture$z1.sims[, 4:6])
sims.culture.tau <- rbind(mediate.culture$tau.sims[, 4:6])


# Cultural threat ACME
# As noted in Online Appendix G (fn. 16)
mean(apply(sims.culture.acme, 1, sum)) # ACME
quantile(apply(sims.culture.acme, 1, sum), c(0.025, 0.975)) # ACME CIs

# ADE
mean(apply(sims.culture.ade, 1, sum))
quantile(apply(sims.culture.ade, 1, sum), c(0.025, 0.975))

# Total effect
mean(apply(sims.culture.tau, 1, sum))
quantile(apply(sims.culture.tau, 1, sum), c(0.025, 0.975))

# Sensitivity analysis
model.culture.outcome.ols <- lm(allowRefugees ~ afterParis + cultureIndex +
                                                gender + ageGroup + education + motherTongue + region + selfPlacement,
                                data = subset(D, allSocio == 1 & !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.culture.outcome.ols)

# Mediation
mediate.culture.ols <- mediate(model.culture, model.culture.outcome.ols, boot = FALSE, treat = "afterParis",
                               mediator = "cultureIndex", sims = n_sims)

mediate.culture.ols.sens <- medsens(mediate.culture.ols, rho.by = 0.01,
                                    effect.type = "indirect")
# Correlation (with unobserved confounder) at which the AMCE is zero
# (as noted in Online Appendix G)
summary(mediate.culture.ols.sens) # "Rho at which ACME = 0: -0.75"

# Figure A3 Panel Bottom Right
# Note: Panels in A3 are assembled in Illustrator
pdf("Graphs/Figure_A3_Bottom_Right-Raw.pdf", 3.5, 3.9, useDingbats = FALSE) # Shrink to 3.25 height
plot(mediate.culture.ols.sens, main = "Culture",
     cex.main = 0.9, cex.axis = 0.85)
dev.off()


# Multi-mediator setup (Table A12 ACME)
model.culture.multimed <- multimed(outcome = "allowRefugeesBinary",
                                   med.main = "cultureIndex",
                                   med.alt = "securityIndex",
                                   treat = "afterParis",
                                   covariates = c("gender", "ageGroup", "education", "motherTongue", "region", "selfPlacement"),
                                   data = subset(D, allSocio == 1 & !is.na(cultureIndex) & !is.na(securityIndex) &
                                                    !is.na(allowRefugeesBinary) & period %in% c(0, 1)),
                                   sims = n_sims)

# Table A12 (Cultural threat ACME)
# As noted in Online Appendix G
summary(model.culture.multimed) # ACME ("ACME (average)")
model.culture.multimed$d.ave.ci[, 1] # Table A12 Cultural threat ACME CIs

# Multi-mediator sensitivity analysis (Not shown)
model.culture.multimed.sens <- multimed(outcome = "allowRefugees",
                                        med.main = "cultureIndex",
                                        med.alt = "securityIndex",
                                        treat = "afterParis",
                                        covariates = c("gender", "ageGroup", "education", "motherTongue", "region", "selfPlacement"),
                                        conf.level = 0.90,
                                        data = subset(D, allSocio == 1 & !is.na(cultureIndex) & !is.na(securityIndex) &
                                                         !is.na(allowRefugeesBinary) & period %in% c(0, 1)),
                                        sims = n_sims)
summary(model.culture.multimed.sens) # Not shown

# Not shown
plot(model.culture.multimed.sens,
     xlim = c(0, 0.3), cex.main = 0.9, cex.axis = 0.85,
     type = c("R2-total"), tgroup = "average", effect.type = "indirect")



##### Economic threat causal mediation regressions (not shown)

# Mediator regression
model.economy <- lm(economyIndex ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                    data = subset(D, !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
# Outcome regression
model.economy.outcome <- polr(allowRefugeesFactor ~ afterParis + economyIndex +
                                                     gender + ageGroup + education + motherTongue + region + selfPlacement,
                              data = subset(D, !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.economy.outcome)

# Mediation
mediate.economy <- mediate(model.economy, model.economy.outcome, boot = TRUE, treat = "afterParis",
                           mediator = "economyIndex", sims = n_sims)
# saveRDS(mediate.economy, "Models/mediate.economy.rds")
# mediate.economy <- readRDS("Models/mediate.economy.rds")



# Simulation results (4:6 represent Slightly/Somewhat/Strongly agree)
# Averaging over z = 0 and z =1
sims.economy.acme <- rbind(mediate.economy$d0.sims[, 4:6],
                           mediate.economy$d1.sims[, 4:6])
sims.economy.ade <- rbind(mediate.economy$z0.sims[, 4:6],
                          mediate.economy$z1.sims[, 4:6])
sims.economy.tau <- rbind(mediate.economy$tau.sims[, 4:6])

# ACME
mean(apply(sims.economy.acme, 1, sum))
quantile(apply(sims.economy.acme, 1, sum), c(0.025, 0.975))

# ADE
mean(apply(sims.economy.ade, 1, sum))
quantile(apply(sims.economy.ade, 1, sum), c(0.025, 0.975))

# Total effect
mean(apply(sims.economy.tau, 1, sum))
quantile(apply(sims.economy.tau, 1, sum), c(0.025, 0.975))

# Sensitivity analysis
model.economy.outcome.ols <- lm(allowRefugees ~ afterParis + economyIndex +
                                                 gender + ageGroup + education + motherTongue + region + selfPlacement,
                                data = subset(D, !is.na(allowRefugeesFactor) & period %in% c(0, 1)))
summary(model.economy.outcome.ols)
# Mediation
mediate.economy.ols <- mediate(model.economy, model.economy.outcome.ols, boot = FALSE, treat = "afterParis",
                               mediator = "economyIndex", sims = n_sims)

mediate.economy.ols.sens <- medsens(mediate.economy.ols, rho.by = 0.01,
                                    effect.type = "indirect")
summary(mediate.economy.ols.sens)

# Not shown
plot(mediate.economy.ols.sens, main = "Economy",
     cex.main = 0.9, cex.axis = 0.85)


# Table A11
# Outcome regression models for mediation analysis
stargazer(model.anxiety.outcome,
          model.sympathy.outcome,
          model.security.outcome,
          model.culture.outcome,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = c("Anxiety", "Sympathy", "Security", "Culture"),
          model.names = FALSE,
          ord.intercepts = TRUE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Anxiety", 
                               "Sympathy", 
                               "Security", 
                               "Culture", 
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "\\midrule$\\kappa_1$",
                               "$\\kappa_2$",
                               "$\\kappa_3$",
                               "$\\kappa_4$",
                               "$\\kappa_5$"))



##########
########## POLITICAL MOBILIZATION GRAPHS AND REGRESSION RESULTS
##########

# Figure A2
pdf("Graphs/Figure_A2-Raw.pdf", 3.25, 2.5, useDingbats = FALSE)
ggplot(subset(G, variable %in% c("Willingness to contact MP\nregarding refugee resettlement")), aes(x = day, y = value, size = treatment, colour = treatment, fill = treatment)) +
  my.theme(bordersize = 1, base_size = 10) +
  # facet_wrap(~ variable, scales = "free", ncol = 1) +
  labs(x = "Days from attack", y = "Proportion would contact MP", fill = "", colour = "") +
  coord_cartesian(xlim = c(-2.55, 20.10), ylim = c(0.37, 0.505), expand = FALSE) +
  scale_x_continuous(limits = c(-2.55, 20.10), expand = c(0, 0),
                     breaks = c(0, 5, 10, 15, 20, 25, 30), labels = c("Paris\nattacks", "5", "10", "15", "20", "25", "30")) +
  scale_y_continuous(breaks = seq(0.3, 0.5, by = 0.04)) +
  geom_vline(xintercept = 0, size = 0.25, linetype = 2, colour = "black") +
  geom_smooth(size = 3, fullrange = FALSE, colour = "grey90", method = "lm", formula = y ~ x + I(x^2), se = FALSE) +
  geom_linerange(aes(ymax = value + contactMPSE * qnorm(0.95), ymin = value - contactMPSE * qnorm(0.95)), size = 0.3) +
  geom_point(shape = 21, size = 2.5) +
  scale_size_manual(values = c("Before attack" = 2.5, "After attack" = 2.25)) +
  scale_fill_manual(values = c("Before attack" = "white", "After attack" = "black")) +
  scale_colour_manual(values = c("Before attack" = "black", "After attack" = "black")) +
  theme(legend.position = "none",
        plot.margin = unit(c(0.5, 0.35, 0.25, 0.25), "lines"))
dev.off()


# Table A6 Model (1)
model.contactMP1 <- glm(contactMP ~ afterParis +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue +
                                    region +
                                    selfPlacement
                                    ,
                                    family = binomial,
                                    data = subset(D, period %in% c(0, 1)))
# summary(model.contactMP1)

# p-value of the effect of the attacks on willingness to contact an MP
# As presented in the main article and in Online Appendix D
coef(summary(model.contactMP1))["afterParis", 4]

# p-value for the unadjusted/raw data
# As presented in Online Appendix D
svychisq(~ contactMP + afterParis,
         design = subset(surveyDesign, period %in% c(0, 1)))


# Table A9 Model (2)
model.contactMP1.immediate.responders <- glm(contactMP ~ afterParis +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue +
                                    region +
                                    selfPlacement
                                    ,
                                    family = binomial,
                                    data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
summary(model.contactMP1.immediate.responders)

# # Table A19 Model (3) (interactive)
model.contactMP1.factorial1 <- glm(contactMP ~ afterParis * treatment1 +
                                    afterParis * treatment2 +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue +
                                    region +
                                    selfPlacement,
                                    family = binomial,
                                    data = subset(D, period %in% c(0, 1)))
summary(model.contactMP1.factorial1)

# Table A19 Model (2) (additive)
model.contactMP1.factorial2 <- glm(contactMP ~ afterParis +
                                               treatment1 +
                                               treatment2 +
                                               gender +
                                               ageGroup +
                                               education +
                                               motherTongue +
                                               region +
                                               selfPlacement,
                                    family = binomial,
                                    data = subset(D, period %in% c(0, 1)))
summary(model.contactMP1.factorial2)


# Table A6 Model (2)
model.contactMP2 <- glm(contactMP ~ I(allowRefugees <= 3) +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue + 
                                    region +
                                    selfPlacement,
                                    family = binomial,
                                    data = subset(D, period == 0))
summary(model.contactMP2)


# p-value of the differences between supporters and opponents of refugee
# resettlement before the attacks
# (as presented in the main article)
coef(summary(model.contactMP2))["I(allowRefugees <= 3)TRUE", 4]


# Table A9 Model (4)
model.contactMP2.immediate.responders <- glm(contactMP ~ I(allowRefugees <= 3) +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue + 
                                    region +
                                    selfPlacement,
                                    family = binomial,
                                    data = subset(D, immediate_responder == 1 & period == 0))
summary(model.contactMP2.immediate.responders)


# Table A6 Model (3)
model.contactMP3 <- glm(contactMP ~ I(allowRefugees <= 3) +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue + 
                                    region +
                                    selfPlacement,
                                    family = binomial,
                                    data = subset(D, period == 1))
summary(model.contactMP3)

# p-value of the differences between supporters and opponents of refugee
# resettlement after the attacks
# (as presented in the main article)
coef(summary(model.contactMP3))["I(allowRefugees <= 3)TRUE", 4]


# Table A9 Model (6)
model.contactMP3.immediate.responders <- glm(contactMP ~ I(allowRefugees <= 3) +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue + 
                                    region +
                                    selfPlacement,
                                    family = binomial,
                                    data = subset(D, immediate_responder == 1 & period == 1))
summary(model.contactMP3.immediate.responders)


# Table A6 Model (4)
model.contactMP4 <- glm(contactMP ~ I(allowRefugees <= 3) * afterParis +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue +
                                    region +
                                    selfPlacement
                                    ,
                                    family = binomial,
                                    data = subset(D, period %in% c(0, 1)))
summary(model.contactMP4)

# p-value of the difference in differences of willingness to contact an MP
# about refugee resettlement between supporters and opponents
# (as presented in the main article)
coef(summary(model.contactMP4))["I(allowRefugees <= 3)TRUE:afterParis", 4]

# Table A9 Model (8)
model.contactMP4.immediate.responders <- glm(contactMP ~ I(allowRefugees <= 3) * afterParis +
                                    gender +
                                    ageGroup +
                                    education +
                                    motherTongue +
                                    region +
                                    selfPlacement
                                    ,
                                    family = binomial,
                                    data = subset(D, immediate_responder == 1 & period %in% c(0, 1)))
summary(model.contactMP4.immediate.responders)


# Sequential-G estimation for estimation of the effect on respondents'
# willingness to contact an MP about refugee resettlement

# As noted in the aticle, there are two separate branches of the survey such
# that respondents received either questions regarding threat, or those
# regarding emotions. To allow both to be accounted for, we first multiply
# impute the missing values in each set of respondents, after which we apply
# sequential-g estimation to estimate the effect of the attacks on willingness
# to contact an MP among supporters and opponents of resettlement
MI_raw <- D[D$period %in% 0:1, c("allowRefugees", "contactMP", "donateIntegrate",
                                 "weight", "period", "afterParis",
                                 "sympathyIndex", "anxietyIndex",
                                 "securityIndex", "cultureIndex", "economyIndex",
                                 "gender", "ageGroup", "education",
                                 "motherTongue", "region", "selfPlacement")]

set.seed(1)
MI_out <- amelia(MI_raw, m = 1,
                 idvars = c("weight", "period"),
                 ords = c("allowRefugees", "selfPlacement"),
                 nom = c("contactMP", "donateIntegrate", "afterParis",
                         "gender", "ageGroup", "education", "motherTongue",
                         "region"))

MI <- MI_out$imputations[[1]]

MI$allowRefugeesBinary <- NA
MI$allowRefugeesBinary[MI$allowRefugees >= 4] <- 1
MI$allowRefugeesBinary[MI$allowRefugees <= 3] <- 0

MI$opposeRefugeesBinary <- NA
MI$opposeRefugeesBinary[MI$allowRefugees >= 4] <- 0
MI$opposeRefugeesBinary[MI$allowRefugees <= 3] <- 1

# For those who oppose refugees (i.e. allowRefugeesBinary == 0)
# The coefficient on afterParis is the effect among opponents as noted in
# the main article
model.contactMPOppose <- sequential_g(contactMP ~ afterParis + gender + ageGroup + education + region + motherTongue + selfPlacement |
                                      anxietyIndex + sympathyIndex + securityIndex + cultureIndex + economyIndex |
                                      allowRefugeesBinary + I(allowRefugeesBinary*afterParis), data = subset(MI, period %in% 0:1))
summary(model.contactMPOppose)
# Effect of the attacks on willingess to contact an MP among opponents
# (as noted in the main article regarding sequential-g estimates)
coef(summary(model.contactMPOppose))["afterParis", 1]*100
# Confidence intervals
round(confint(model.contactMPOppose)[2, ], 3)*100

# For those who support refugees (i.e. opposeRefugeesBinary == 0)
# The coefficient on afterParis is the effect among opponents as noted in
# the main article
model.contactMPSupport <- sequential_g(contactMP ~ afterParis + gender + ageGroup + education + region + motherTongue + selfPlacement |
                                      anxietyIndex + sympathyIndex + securityIndex + cultureIndex + economyIndex |
                                      allowRefugeesBinary + I(opposeRefugeesBinary*afterParis), data = subset(MI, period %in% 0:1))
summary(model.contactMPSupport)
# Effect of the attacks on willingess to contact an MP among supporters
# (as noted in the main article regarding sequential-g estimates)
coef(summary(model.contactMPSupport))["afterParis", 1]*100
# Confidence intervals
round(confint(model.contactMPSupport)[2, ], 3)*100


# Difference in differences in willingess to contact an MP among supporters and
# opponents of refugees resettlement
# (as noted in the main article regarding sequential-g estimation)
(coef(summary(model.contactMPOppose))["afterParis", 1] -
 coef(summary(model.contactMPSupport))["afterParis", 1]) * 100


# Raw data estimate among opponents (not shown)
svymean(~ contactMP, design = subset(surveyDesign, allowRefugees <= 3 & period %in% c(0)), na.rm = TRUE) -
svymean(~ contactMP, design = subset(surveyDesign, allowRefugees <= 3 & period %in% c(1)), na.rm = TRUE)

# Raw data estimate among supporters (not shown)
svymean(~ contactMP, design = subset(surveyDesign, allowRefugees >= 4 & period %in% c(0)), na.rm = TRUE) -
svymean(~ contactMP, design = subset(surveyDesign, allowRefugees >= 4 & period %in% c(1)), na.rm = TRUE)

# Table A6
stargazer(model.contactMP1, model.contactMP2,
          model.contactMP3, model.contactMP4,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = TRUE,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Oppose refugees", 
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Oppose refugees $\\times$ Paris attacks",
                               "Intercept"))

# Table A9
stargazer(model.contactMP1,
          model.contactMP1.immediate.responders,
          model.contactMP2,
          model.contactMP2.immediate.responders,
          model.contactMP3,
          model.contactMP3.immediate.responders,
          model.contactMP4,
          model.contactMP4.immediate.responders,
          title = "",
          style = "apsr",
          column.sep.width = "0pt",
          star.cutoffs = c(0.05, 0.01, 0.001),
          dep.var.labels.include = TRUE,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c("Paris attacks",
                               "Oppose refugees", 
                               "Female",
                               "Age 30-39",
                               "Age 40-49",
                               "Age 50-64",
                               "Age 65+",
                               "College",
                               "University degree",
                               "Francophone",
                               "Other language",
                               "Ontario",
                               "Quebec",
                               "West",
                               "Political ideology",
                               "Oppose refugees $\\times$ Paris attacks",
                               "Intercept"))


# Graph the predicted probabilities of willingness to contact an MP before
# and after the attacks (for Figure 4)
get_sim_estimates <- function(fitted_model, n_sims, allow_refugees) {

  # Simulate coefficients
  set.seed(1)
  coefs_sim <- rmvnorm(n_sims, coef(fitted_model), vcov(fitted_model))

  # Model matrices with predictors for the sensitive-item and misreport sub-models
  data_sim <- subset(D, !is.na(weight))[-which(names(D) == "allowRefugees"), ]
  data_sim$allowRefugees <- allow_refugees
  model_matrix <- data.frame(model.matrix(fitted_model, data = data_sim))
  weight_sim <- data_sim$weight[match(rownames(model_matrix), rownames(data_sim))]

  sims <- sapply(1:n_sims, function(x) weighted.mean(plogis(coefs_sim[x, ] %*% t(model_matrix)), weight_sim))
  return(data.frame(value = mean(sims),
                    lower = quantile(sims, 0.05),
                    upper = quantile(sims, 0.95)))

}

G3 <- data.frame(group = rep(c("Support resettlement", "Oppose resettlement"), each = 2),
                 treatment = rep(c("Before attacks", "After attacks"), 2),
                 value = NA, lower = NA, upper = NA)

# Simulate predicted probabilities (1 minute to run)
set.seed(1)
G3[, c("value", "lower", "upper")] <- rbind(get_sim_estimates(model.contactMP2, 5000, 6)*100,
                                            get_sim_estimates(model.contactMP3, 5000, 6)*100,
                                            get_sim_estimates(model.contactMP2, 5000, 1)*100,
                                            get_sim_estimates(model.contactMP3, 5000, 1)*100)

G3$treatment <- factor(G3$treatment, levels = c("Before attacks", "After attacks"))
G3$group <- factor(G3$group, levels = c("Support resettlement", "Oppose resettlement"))
G3$group2 <- as.character(G3$group)
G3$group2[G3$group2 == "Support resettlement"] <- "Support\nresettlement"
G3$group2[G3$group2 == "Oppose resettlement"] <- "Oppose\nresettlement"
G3$group2 <- factor(G3$group2, levels = c("Support\nresettlement", "Oppose\nresettlement"))

# Difference in predicted probabilities of expressing willingness to contact
# an MP about refugee resettlement before the attacks among opponents and
# supporters of resettlement
# (as noted in the main article, 3 percentage points)
G3$value[G3$group == "Support resettlement" & G3$treatment == "Before attacks"] -
G3$value[G3$group == "Oppose resettlement" & G3$treatment == "Before attacks"]

# Difference in predicted probabilities of expressing willingness to contact
# an MP about refugee resettlement before and after the attacks among opponents
# and supporters of resettlement
# (as noted in the main article, 11 versus 2 percentage points)
G3$value[G3$group == "Support resettlement" & G3$treatment == "Before attacks"] -
G3$value[G3$group == "Support resettlement" & G3$treatment == "After attacks"]

G3$value[G3$group == "Oppose resettlement" & G3$treatment == "Before attacks"] -
G3$value[G3$group == "Oppose resettlement" & G3$treatment == "After attacks"]



# Figure 4
pdf("Graphs/Figure_4-Raw.pdf", 3.5, 2.3, useDingbats = FALSE)
ggplot(G3, aes(x = treatment, y = value,
       colour = group, fill = group)) +
  my.theme(borderless = 2, base_size = 9) +
  labs(color = "", fill = "", size = "", y = "Pr(Contacting MP)", x = "") +
  coord_cartesian(ylim = c(0, 65), expand = FALSE) +
  scale_y_continuous(breaks = c()) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9), width = 0.7, size = 0.4) +
  geom_linerange(aes(ymax = upper, ymin = lower), size = 0.5,
                 position = position_dodge(width = 0.9), color = "black") +
  geom_hline(yintercept = 0, size = 1) +
  geom_text(aes(y = upper + 6, label = round(value)), position = position_dodge(width = 0.9), colour = c("black", "black", "black", "black"), size = 3.2) +
  scale_color_manual(values = c("Support resettlement" = "black",
                                "Oppose resettlement" = NA)) +
  scale_fill_manual(values = c("Support resettlement" = "white",
                               "Oppose resettlement" = "grey80")) +
  theme(legend.position = "top",
        legend.direction = "vertical",
        plot.margin = unit(c(-1.25, 4.5, -0.5, 1), "lines"))
dev.off()



# Online Appendix C
# Bon Ferroni corrected p-values (multiple comparisons correction)
p_all <- c(anxiety = summary(model.anxiety)$coefficients["afterParis", 4],
           security = summary(model.security)$coefficients["afterParis", 4],
           culture = summary(model.culture)$coefficients["afterParis", 4],
           sympathy = summary(model.sympathy)$coefficients["afterParis", 4],
           economy = summary(model.economy)$coefficients["afterParis", 4],
           support_refugees = coeftest(model.refugees1)["afterParis", 4],
           contact_mp = summary(model.contactMP1)$coefficients["afterParis", 4],
           contact_mp_oppose_refugees = summary(model.contactMPOppose)$coefficients["afterParis", 4])

bon_ferroni <- p.adjust(p_all, method = "bonferroni")

# Summary of bon ferroni-corrected p-values
# (as noted in Online Appendix C "Multiple comparisons corrections")
data.frame(unadjusted = ifelse(p_all < 0.001, "p < 0.001",
                               ifelse(p_all < 0.001, "p < 0.01",
                                      ifelse(p_all < 0.05, "p < 0.05", "Not sig."))),
           adjusted = ifelse(bon_ferroni < 0.001, "p < 0.001",
                             ifelse(bon_ferroni < 0.001, "p < 0.01",
                                    ifelse(bon_ferroni < 0.05, "p < 0.05", "Not sig."))))


# Online Appendix I

# Likelihood ratio tests to check whether the religion and location conditions
# jointly modify the effect of the Paris attacks on each main
# outcome of interest. In no case does including these moderators improve
# model fit

# Anxiety likelhood ratio test
# (Table A14 Model (2) versus Model (3))
# As noted in Online Appendix I
lrtest(model.anxiety.factorial1,
       model.anxiety.factorial2)$`Pr(>Chisq)`[2]

# Sympathy likelhood ratio test
# (Table A15 Model (2) versus Model (3))
# As noted in Online Appendix I
lrtest(model.sympathy.factorial1,
       model.sympathy.factorial2)$`Pr(>Chisq)`[2]

# Security threat likelhood ratio test
# (Table A16 Model (2) versus Model (3))
# As noted in Online Appendix I
lrtest(model.security.factorial1,
       model.security.factorial2)$`Pr(>Chisq)`[2]

# Cultural threat likelhood ratio test
# (Table A17 Model (2) versus Model (3))
# As noted in Online Appendix I
lrtest(model.culture.factorial1,
       model.culture.factorial2)$`Pr(>Chisq)`[2]

# Support for resettlement likelhood ratio test
# (Table A18 Model (2) versus Model (3))
# As noted in Online Appendix I
lrtest(model.refugees1.factorial1,
       model.refugees1.factorial2)$`Pr(>Chisq)`[2]

# Willing to contact an MP likelhood ratio test
# (Table A19 Model (2) versus Model (3))
# As noted in Online Appendix I
lrtest(model.contactMP1.factorial1,
       model.contactMP1.factorial2)$`Pr(>Chisq)`[2]

covariate_labels <- c("Paris attacks",
                     "Condition 1: Christian",
                     "Condition 1: Muslim",
                     "Condition 2: Community",
                     "Female",
                     "Age 30-39",
                     "Age 40-49",
                     "Age 50-64",
                     "Age 65+",
                     "College",
                     "University degree",
                     "Francophone",
                     "Other language",
                     "Ontario",
                     "Quebec",
                     "West",
                     "Political ideology",
                     "Condition 1: Christian $\\times$ Paris attacks",
                     "Condition 1: Muslim $\\times$ Paris attacks",
                     "Condition 2: Community $\\times$ Paris attacks")

# Table A14
stargazer(model.anxiety, model.anxiety.factorial2, model.anxiety.factorial1,
          title = "", style = "apsr", column.sep.width = "0pt",
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = NULL,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c(covariate_labels, "Constant"))

# Table A15
stargazer(model.sympathy, model.sympathy.factorial2, model.sympathy.factorial1,
          title = "", style = "apsr", column.sep.width = "0pt",
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = NULL,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c(covariate_labels, "Constant"))

# Table A16
stargazer(model.security, model.security.factorial2, model.security.factorial1,
          title = "", style = "apsr", column.sep.width = "0pt",
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = NULL,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c(covariate_labels, "Constant"))

# Table A17
stargazer(model.culture, model.culture.factorial2, model.culture.factorial1,
          title = "", style = "apsr", column.sep.width = "0pt",
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = NULL,
          model.names = FALSE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c(covariate_labels, "Constant"))

# Table A18
stargazer(model.refugees1, model.refugees1.factorial2, model.refugees1.factorial1,
          title = "", style = "apsr", column.sep.width = "0pt",
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = NULL,
          model.names = FALSE,
          ord.intercepts = TRUE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c(covariate_labels,
                               "\\midrule$\\kappa_1$",
                               "$\\kappa_2$",
                               "$\\kappa_3$",
                               "$\\kappa_4$",
                               "$\\kappa_5$"))

# Table A19
stargazer(model.contactMP1, model.contactMP1.factorial2, model.contactMP1.factorial1,
          title = "", style = "apsr", column.sep.width = "0pt",
          dep.var.labels.include = FALSE,
          star.cutoffs = c(0.05, 0.01, 0.001),
          column.labels = NULL,
          model.names = FALSE,
          ord.intercepts = TRUE,
          omit.stat = c("ll", "aic"),
          covariate.labels = c(covariate_labels,
                               "Constant"))





# Online Appendix H

# Set up the max and min y-axis limits for Figure 2 (bit of a hack for ggplot using geom_blank())
G$max <- G$min <- as.numeric(NA)
G$min[G$variable == "Retrospective evaluation\nof personal financial situation"] <- -0.3
G$max[G$variable == "Retrospective evaluation\nof personal financial situation"] <- 0.3

G$min[G$variable == "Retrospective evaluation\nof national economy"] <- -0.3
G$max[G$variable == "Retrospective evaluation\nof national economy"] <- 0.3


# Figure A5
pdf("Graphs/Figure_A5-Raw.pdf", 6, 2.6, useDingbats = FALSE)
ggplot(subset(G, variable %in% c("Retrospective evaluation\nof personal financial situation",
                                 "Retrospective evaluation\nof national economy")),
       aes(x = day, y = value, size = treatment, colour = treatment, fill = treatment)) +
  my.theme(bordersize = 1, base_size = 9) +
  facet_wrap(~ variable, scales = "free", ncol = 2) +
  labs(x = "Days from attack", y = "Economic evaluation (std. units)", fill = "", colour = "") +
  # coord_cartesian(xlim = c(-2.55, 20.10), expand = FALSE) +
  scale_x_continuous(limits = c(-2.55, 20.10), expand = c(0, 0), breaks = c(0, 5, 10, 15, 20, 25, 30)) +
  scale_y_continuous(limits = c(-0.304, 0.304), expand = c(0, 0), breaks = seq(-0.3, 0.3, by = 0.1), labels = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3)) +
  geom_vline(xintercept = 0, size = 0.25, linetype = 2, colour = "black") +
  geom_smooth(size = 3.5, fullrange = FALSE, colour = "grey90", method = "lm", formula = y ~ x + I(x^2), se = FALSE) +
  geom_linerange(data = subset(G, variable == "Retrospective evaluation\nof personal financial situation"), aes(ymax = value + economyRetrospectiveSelfSE * qnorm(0.95), ymin = value - economyRetrospectiveSelfSE * qnorm(0.95)), size = 0.3) +
  geom_linerange(data = subset(G, variable == "Retrospective evaluation\nof national economy"), aes(ymax = value + economyRetrospectiveCanadaSE * qnorm(0.95), ymin = value - economyRetrospectiveCanadaSE * qnorm(0.95)), size = 0.3) +
  geom_point(shape = 21, size = 2.5) +
  geom_blank(aes(y = max)) +
  geom_blank(aes(y = min)) +
  scale_size_manual(values = c("Before attack" = 2.5, "After attack" = 2.25)) +
  scale_fill_manual(values = c("Before attack" = "white", "After attack" = "black")) +
  scale_colour_manual(values = c("Before attack" = "black", "After attack" = "black")) +
  theme(legend.position = "none",
        plot.margin = unit(c(0.5, 0.35, 0.25, 0.25), "lines"))
dev.off()

# Effect of the attacks on retrospective economic evaluations in Canada
# As noted in Online Appendix H
model.retroCanada <- glm(economyRetrospectiveCanada ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                         data = subset(D, period %in% c(0, 1)))
# summary(model.retroCanada)
coef(summary(model.retroCanada))["afterParis", 4] # 0.82

# Effect of the attacks on retrospective economic evaluations for oneself
# As noted in Online Appendix H
model.retroSelf <- glm(economyRetrospectiveSelf ~ afterParis + gender + ageGroup + education + motherTongue + region + selfPlacement,
                       data = subset(D, period %in% c(0, 1)))
# summary(model.retroSelf)
coef(summary(model.retroSelf))["afterParis", 4] # 0.16


